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Abstract. The resonance states and the decay dynamics of the nonhnear Schrodinger 
(or Gross-Pitaevskii) equation are studied for a simple, however flexible model system, 
the double delta-shell potential. This model allows analytical solutions and provides 
insight into the influence of the nonlinearity on the decay dynamics. The bifurcation 
scenario of the resonance states is discussed, as well as their dynamical stability 
properties. A discrete approximation using a biorthogonal basis is suggested which 
allows an accurate description even for only two basis states in terms of a nonlinear, 
nonhcrmitian matrix problem. 

PACS numbers: 03.65Ge, 03.65Nk, 03.75-b 
1. Introduction 

In quantum dynamics metastable states can be conveniently described as resonance 
states, i.e. eigenstates of the Schrodinger equation with complex eigenvalues. Such 
resonances can be efficiently calculated by complex scaling methods [1] or matrix 
truncation techniques for periodic lattices [2] and found numerous applications. 
The recent progress of the physics of Bose-Einstein condensates (BEC) stimulated 
investigations of the role of resonances for such systems, as for instance escape from 
a potential well or decay by loss of condensate particles as it was realized in a recent 
experiment with ultracold molecules in a one-dimensional lattice [3]. Alternative 
implementations of open systems with interactions include experiments with optical 
wave guide arrays [4-6]. 

Here we will consider a (quasi) one-dimensional configuration described by the 
nonlinear Schrodinger equation (NLSE) or Gross-Pitaevskii equation (GPE) for the 
macroscopic condensate wavefunction 

where the nonlinear parameter g describes the self-interaction. First applications of 
complex energy resonance states extending the complex scaling method to the nonlinear 
case investigated the decay of a condensate trapped in the celebrated model potential 
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V{x) ~ (x^ — /3)e~"^ [7^9]- However, the definition of a resonance is somewliat 
ambiguous in tlie nonlinear case; alternative descriptions based on amplitude ratios 
in the inner and outer potential region have also been proposed and applied to simple 
model systems allowing an analytical treatment [10]. One of these models is the delta- 
shell potential discussed in section O 

Related studies explore resonance phenomena observed in transport problems of 
BECs, for example the transmission though a potential barrier or across a potential 
well [11-14] which can be explained in terms of an underlying resonance state [13,14]. 
Complex resonances have also been used to describe a BEC in accelerated optical 
lattices by means of nonlinear Wannier-Stark states [15, 16]. Here the potential in 
([1]) is of the form Vp{x) + Fx, where Vp{x) is periodic in space. Such systems can 
also be described in terms of Wannier functions localized on the potential minima, the 
'sites'. In a single-band approximation this leads to a discrete nonlinear Schrodinger (or 
Gross-Pitaevskii) equation, well known as the discrete self-trapping equation (DST) [17], 
which is even of interest in its most simple form of only two sites. Alternatively, 
this equation can be derived starting from a many-particle description of a BEC on 
a discrete lattice by a Bose-Hubbard Hamiltonian, which leads again to the DST 
equation in the mean-field limit. In such a description, however, the decay has been 
neglected. In can be re-introduced again by opening the system. This can be done in 
various ways, purely phenomenologically by introducing complex site energies describing 
decay [18,19] or more sophisticatedly by taking explicitly the coupling to an environment 
into account. Recent studies comparing the full many particle dynamics with the mean- 
field approximation consider the two site system (the open 'dimer') phenomenologically 
[18, 19] or using a master equation for the coupling to the environment [20-22] (see 
also [23-25] for a study of the related Lipkin-Meshkov-Glick model). It should be 
noted that even the resulting nonhermitian nonlinear two-level system shows an intricate 
crossing scenario as discussed, e.g. , in [26]. Two-level systems are often used to model 
double-well potentials realized in various BEC experiments [27-31]. 

The present paper is devoted to a detailed analysis of a simple model system, 
the decay behavior of a BEC in a double delta-shell potential. This open double- 
well system is on the one hand simple enough to allow an analytic treatment and 
closed form approximations and on the other hand it is flexible enough to investigate 
characteristic phenomena observed in nonlinear double-well dynamics, as self-trapping 
and the appearance of new eigenvalues through a saddle node bifurcation [32] and their 
modification due to the decay. 

The paper is organized as follows: In section [2] we first discuss the nonlinear 
single delta-shell potential and derive simple analytic approximations for the resonance 
position and decay rate which are compared with exact numerical results. These 
techniques are in section [3] extended to the double delta-shell potential where the 
bifurcation scenario of nonlinear resonance states is analyzed as well as the decay 
dynamics. A discrete basis set expansion is used in section 13.31 to reduce the system to 
a finite nonlinear, nonhermitian matrix problem, which yields reasonable results even 
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for a two state approximation. A Bogoliubov-de-Gennes analysis in section [231 provides 
information about the stability of the resonance states. Additional material concerning 
computational details is present in an appendix. 



2. Single delta-shell 

At first we consider the case of an open single well, namely the so-called delta-shell 
potential 

V{x) = [ ''-^ (2) 

^ ' \ {hym)X6{x-a) x>0 ^ ' 

with A > 0, a > and repulsive interaction g > 0. This potential and its generalization 
to three dimensions have been investigated in detail in the context of the linear 
Schrodinger equation [33-36]. In the context of the GPE it has been considered in [10] 
in which resonance positions and widths are extracted from real valued wave functions 
ip{x) and a resonance is characterized by a maximum of the probability density 
inside the potential region < x < a. 

Here we consider nonlinear resonance states, i.e. eigenstates of the time-independent 
NLSE with complex eigenvalues /i — ir/2 

+ it^- ir/2 - V)ij - = (3) 

and purely outgoing (Siegert) boundary conditions (for details see [14]) which enables 
us to determine both position fi and decay width r/2 of the resonances and to 
construct an approximation that allows an analytical treatment. Note that the states 
ip{x, t) = exp [— i(yU — ir/2)t/h] tplx) do not satisfy the time-dependent NLSE since their 
norm is not constant. Instead they provide an adiabatic approximation to the actual 
time-dependence (see [9,37] and section [3l2l) . 

The NLSE ([3]) with Siegert boundary conditions is solved by the ansatz 

, , , J Isn{gx\p) < X < a 

^ 0^ = S ^ ikx (4) 

1 C e'"^ X > a ^ ^ 

with a Jacobi elliptic sn-function inside and an outgoing plane wave outside the potential 
well. The parameters have to satisfy < p < 1, 



and 

gm 

We fix the norm of the wavefunction to unity inside the potential region < x < a: 



(6) 



1 = r |^(x)|Mx = ^{ga- E{ga\p)) , (7) 
Jo gm 
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where E{u\p) is the incomplete eUiptic integral of the second kind. The matching 
conditions at x = a, ip{a+) = ip{a—) = ip{a), ip'{a+) — ip'{a—) = 2\ip{a), yield 

/sn(^a|p) = Ce''^'^ (8) 
/^cn(^a|p)dn(^a|p) + 2A/sn(^a|p) = ikC e''"' . (9) 

For a strongly repulsive delta-function at x = a we have \k\ <^ A and sn{Qa\p) ~ so 
that the decay of the lowest resonances is weak. Thus we neglect the imaginary part 
— r/2 of the eigenvalue in equation ([5]) and ([9]). For A — oo the wave number g is given 
hj g = 2nK{p)/a {K{p) is the complete elliptic integral of the first kind) so that for 
weak decay we can assume 
2nK{p) 



S 



with 5<0, |(5-a|<^l. Expanding the real part of equation 
we obtain 



(10) 

up to second order in S 
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2Aa + 1 



2Aa 



2nK{p){l+p)a \j \2nK{p){l + p)a J {1 + p)a'^ 

For a given value of the Jacobi elliptic parameter < p < 1 the real part of the 
eigenvalue is determined by the equations ([5]), ( ITOi) and ( ITTi) . From the normalization 
condition ([7]) the interaction strength is given by 

9 = —{Qa- E{ga\p)) . (12) 
m 

The imaginary part of the eigenvalue can be estimated using the Siegert relation [14,38] 
h^k |^(a)P fi^kgp sn^{ga\p) 



r/2 



2m \ip{x)\'^dx 2m ga — E{ga\p) 



(13) 



> 1 




Figure 1. Resonance wave function for tlie ground and first excited state of a delta- 
shell potential (parameters A = 10, a = 1 and nonlinearity g ~ 1; units h ~ m = 1). 
Left: ground state, n = 1, right: first excited state, n = 2. The exact numerical 
solutions obtained by the CAP method (dashed blue) are in excellent agreement with 
the analytical approximations (solid red). 
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Exact results can be obtained from a numerical calculation combining a complex 
absorbing potential (CAP) with a grid relaxation. In these computations we use the 
approximate resonances derived above as starting values. Complex absorbing potentials 
were shown to be equivalent to exterior complex scaling [1] and successfully applied to 
resonances of the GPE [7]. A detailed discussion of our particular implementation as well 
as its application to nonlinear Wannier-Stark resonances will be presented elsewhere [39]. 

As an example, figure [1] shows the wavefunctions for the ground state and the first 
excited state of the delta-shell potential ([2]) with A = 10, a = 1 and a nonlinearity of 
(7 = 1 (units with h = 1 = m are used throughout this paper) and table [1] lists the 
resonance energies. In addition the analytical approximations derived above are given. 
Good agreement is observed. 



n 




AtCAP 


rA/2 


rcAp/2 


1 


5.8950 


5.9043 


0.0759 


0.0751 


2 


19.4138 


19.5065 


0.4810 


0.4680 



Table 1. Resonance energies /i„ — ir„/2 for the two lowest resonances of the 
deha-shell potential (parameters A = 10, a = 1 and nonlinearity g = \) calculated 
analytically (approximation A) and numerically using a grid relaxation method with 
complex absorbing potentials (CAP). 



3. Double delta-shell 

An open double well allows to study the influence of the decay on characteristic nonlinear 
phenomena like self-trapping. Here we investigate a simple case, namely the double- 
delta-shell potential 

, . _ J +00 X < . . 

^"^^^ \{h^/m)[\bS{x-h) + \aS{x-a)] x>0 ^ ^ 

which consists of an infinitely high wall and two repulsive delta barriers with strength 
Aa > 0, Af, > located at a > 6 > 0. This potential possesses two simpler limits: For 
Aa ^ or Af, ^ we recover the single delta-shell potential discussed in section O 
In the limit A^ — > oo the system is closed yielding a coherent nonlinear tunneling 
oscillation between two (asymmetric) potential wells [32,40]. For finite values of Aa we 
therefore observe nonhermitian generalizations of these simple cases. The considerations 
in the following subsection essentially follow those from the preceding section for the 
single delta-shell potential. We first consider stationary resonance states followed by 
a discussion of the djTiamics and a stability analysis. As in the preceding section we 
concentrate on repulsive nonlinearities (7 > 0. 
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3.1. Stationary States 

The NLSE ([3]) with the potential (fT4l) and Siegert boundary conditions is solved by the 
ansatz 

{/isn(f)ix|]?i) < a; < 6 
I2S\i{q2X + 'd\p2) b < X < a (15) 
C e^^"^ X > a 

with 

gm ' gm 

and 0<p<l. Atx = awe obtain the matching conditions 

J2sn(^?2a + ^?b2) = Ce'^» (18) 

/2^2cn(^2a + ^?b2)dn(^2a + ^\P2) + 2A„/2sn(^2a + ^^^2) = ikC e'^" . (19) 

Due to the repulsive delta-function at x = a we have /c <C and sn(^2a + '&\P2) ~ 0. 
As in the previous section we neglect the decay coefficient r/2 in equation ( IlQll . 
Furthermore we assume Q2a + 'd = 4:nK{p2) + 6 ■ a with 5-a<0, |(5-a|<^l. Expanding 
the real part of equation ( fT9i) up to second order in S we obtain 



*-T7^-JfT7^1 +7^^. (20) 



^2(1+^2)0^ V V^'2(l +P2)aV (1+^2)0^ 

so that the phase shift is given by = 4nK{p2) + 6 ■ a — Q2a. 
The matching conditions at x = 6 read 

Jisn(uilpi) = J2sn(u2|p2) , (21) 
/i^icn(Mi|pi)dn(Mi|pl) + 2A;,/isn(Mi|pi) = /2^i2cn(M2|p2)dn(M2|p2) (22) 

where Ui = gib, U2 = Q2b + 'd. By inserting equation fl2T|) into fl22|) and using the 
relations between the squares of the Jacobi elliptic functions we arrive at 

-, — — TTT = 7 — -, ^^--7777^^ sn(M2li02)cn(M2li02)dn(-U2li02) 



+ ^^^sn(n2|p2)' =: F{p2,^^) (23) 
P2 + J- 

so that the parameter pi is given by 

From equations fl20l) - (1241) all relevant quantities are known in terms of p2 and /i. To 
calculate the eigenvalues we solve equation fl22l) for a given value of /i and thus obtain 
up to four solutions for p2 in [0, 1]. 
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m I 



gi {gib - E{gib\pi)) 

+ g2 {g2{a -b)+ E{g2b + i9\p2) - E{g2a + i^\p2)) 
The decay coefficient is again estimated by the Siegert relation 



r/2 



2m |^/;(x)Pdx 



(25) 



(26) 
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Figure 2. Analytically calculated bifurcation diagrams for the parameters Af, = 10, 
Aq = 20, 6=1 and a = 2. Left panel: Chemical potentials. Right panel: Decay rates. 





-0.5 



Figure 3. Wavefunctions of the two lowest states for Ab = 10, Aq = 20, 6 = 1, a = 2. 
Left panel: ground state for g ~ (solid black) and g = 0.5 (dashed red). Right panel: 
first excited state for g = (solid black) and g = 0.5 (dashed green). 

As an example we consider the lowest eigenstates of the potential with parameters 
Aft = 10, Aa = 20, 6 = 1 and a = 2. Figure [2] shows the real and imaginary parts of 
the eigenvalues in dependence of the nonlinear parameter g. The two lowest eigenvalues 
H of the linear {g = 0) system (dashed green and dashed red) increase almost linearly 
with increasing repulsive interaction strength g. The corresponding wavefunctions are 
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Figure 4. Analytically calculated wavefunctions for Ab = 10, \a = 20, 6 = 1, a = 2 
and (7 = 3. Upper left: autochtonous self-trapping state (AuT), lower left: allochtonous 
self-trapping state (AIT), upper right: allochtonous almost antisymmetric state (A1-), 
lower right: autochtonous almost symmetric state (Au-|-). 



Eigenstate 




rA/2 


l^CKP 


rcAp/2 


autochtonous self-trapping state (AuT) 


8.589 


0.00114 


8.588 


0.00114 


allochtonous self-trapping state (AIT) 


8.301 


0.02453 


8.303 


0.02443 


allochtonous almost antisymmetric state (A1-) 


6.999 


0.01279 






autochtonous almost symmetric state (Au-|-) 


6.014 


0.00902 


6.014 


0.00900 



Table 2. Chemical potential and decay rates for the same states as in figure [HI The 
approximate analytical values (A) are compared with numerically exact ones (CAP) 
calculated by a grid relaxation method with complex absorbing potentials. 



shown in figures [3] and H] for g = 0, 0.5 and 3. For a weak nonlinearity (figure [3]) the 
the ground state is almost symmetric and the first excited state almost antisymmetric. 
These states with linear counterpart are referred to as autochtonous states. At (7 ~ 1 
two new eigenvalues appear through a saddle node bifurcation (cf. [32]) which we will 
henceforth call allochtonous states. One of these two states is dynamically stable (dashed 
dotted black) the other unstable (solid blue) (see subsection 13.21) . After the bifurcation 
the state, which was almost antisymmetric before the bifurcation, (dashed green) more 
and more localizes in the left well with increasing interaction, whereas the newly created 
state (dashed dotted black) localizes in the right well (see upper left and lower left panel 
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Figure 5. Analytical results for Ab = 10, \a oo, 6 = 1, a = 2. Upper left 
panel: Bifurcation diagram. The other panels show eigenstates for (7 = 3. Lower left: 
self-trapping state. Upper right:antisyminetric state. Lower right: symmetric state. 



of figure Hj). These two symmetry breaking states are referred to as autochtonous self- 
trapping states (AuT) and allochtonous self-trapping states (AIT), respectively. The 
two remaining states are referred to as the autochtonous almost symmetric state (Au-|-) 
(dashed red) and allochtonous almost antisymmetric state (A1-) (solid blue) (see upper 
right and lower right panel of figure H]) . 

The imaginary parts of the eigenvalues, the decay widths r/2, shown in the right 
panel of[2]grow with increasing interaction strength g for large values of g. This behavior 
can be understood via the Siegert formula (l26i) which predicts a dependence r/2 oc ^/Ji 
if the shape of the wavefunction does not change much with yU (respectively g) and 
furthermore is proportional to g in this parameter range (cf. left panel of figure [2]). 
For small values of g the decay coefficient of the autochtonous self-trapping state (dashed 
green) decreases rapidly with increasing g because of the increasing localization of the 
wavefunction in the left well where the probability for tunneling out of the barrier region 
< X < a is small. 

For comparison we consider the limit Aa — > oo in which the system becomes both 
symmetric and hermitian. Similar systems have recently been considered in a number 
of papers [32,40-42]. Naturally the analytical results obtained in this limit are exact. 
From equation fl20l) we see that in this case we have 5 = and all eigenvalues are real. In 
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X X 

Figure 6. Comparison between the squared magnitudes of the analytically calculated 
eigenstates (solid red) and the numerically exact solution obtained by the CAP method 
(dashed blue) for the same states as in figure [H i.e. upper left: autochtonous self- 
trapping state (AuT). lower left: allochtonous self-trapping state (AIT), upper right: 
allochtonous almost antisymmetric state (A1-), lower right: autochtonous almost 
symmetric state (Au-|-). On the scale of drawing the results of the two methods are 
almost indistinguishable. 
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Figure 7. Overlap of the two allochtonous states for the parameters Af, = 10, Xa = 20, 
b=l,a = 2. 



the upper left panel of figure [5] the chemical potential fi is plotted in dependence of the 
interaction strength g. Compared to the nonhermitian nonsymmetric case considered 
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before the saddle node bifurcation has changed to a pitchfork bifurcation since the 
eigenvalues of the two self-trapping states now coincide due to the symmetry of the 
system (cf. [32]). One of these states is shown in the lower left panel of figure [5l the 
other one is obtained by mirroring at the axis x = b. Hence the respective eigenvalues, 
indicated by the dashed dotted black curve, are degenerate after the bifurcation. The two 
remaining eigenstates (lower right and upper right panel) are now genuinely symmetric, 
respectively antisymmetric. 

As in the case of the open single well we compare our approximate analytical 
results for the resonance solutions of the open double well with numerically exact 
ones calculated using a grid relaxation method with complex absorbing potentials (see 
section [2]). The results are given in figure [6] and table [21 For the two autochtonous states 
we observe good agreement between the analytical approximation and the numerically 
exact solutions. For the AIT state (dashed dotted black) the grid relaxation only 
converges for high values of g (respectively /i) but wherever it converges there is good 
agreement between the analytical approximation and the numerically exact solutions. 
For the Al- state (solid blue) the grid relaxation does not converge at all. It turns out 
(see sections [32] and [331) that this state is dynamically unstable just like the respective 
state in the asymmetric hermitian double well (see [32]). 

At the bifurcation point the real and imaginary parts of the eigenvalues of the 
allochtonous states both coincide. By computing the overlap of the two respective 
wavefunctions (figure [7]) one can demonstrate that the wavefunctions also coincide at 
the bifurcation point. Thus the bifurcation point is an example of an exceptional point 
(see [18,43,44]). 

3.2. Dynamics 

The eigenstates calculated in the preceding subsection have complex eigenvalues and 
are hence subject to decay. Following [9,37] we call an eigenstate dynamically stable if 
its time evolution follows the stationary adiabatic decay behavior given by 



where g^ff = gN with N = \ip{x)\'^dx denotes the effective nonlinear interaction in 
inside the double well potential [10]. In subsection 13.11 the total norm was kept fixed 
at iV = 1 and g was varied. Now we keep g fixed whereas N decreases due to decay. 
We compute the time evolution of the states shown in figure [H using a Crank-Nicholson 
propagation with a predictor-corrector algorithm and absorbing boundaries [13,45]. 

The results are given in figures [HI- [TTl each of which shows a spatio-temporal contour 
plot of the density \ip{x,t)\'^ (left panel), the relative population of the two wells (upper 
right panel) and the decay of the norm (lower right panel) which is calculated directly 
via time-evolution (dashed blue curve) and compared to the stationary decay behavior 
according to equation ( 1271) (solid red curve) calculated from the stationary eigenvalues 
shown in figure 121 The time-evolution of the autochtonous states (figures M and ITTll 
closely follows the stationary decay behavior. Thus these states are dynamically stable. 




(27) 
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t 

Figure 8. Time evolution of the autochtonous self-trapping state (AuT) from 
figure m Left panel: Spatio-temporal contour plot of the density |'(/'(a;)p. Upper 
right panel: relative population of the two wells. iVi and iV2 denote the total norm 
inside the left respectively right well. Lower right panel: Decay of the norm N of the 
wavefunction inside the double barrier calculated via time evolution (dashed blue) and 
using stationary states (solid red). Nq is the norm at t = 0. 
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Figure 9. The same as figure [51 however, for the allochtonous almost antisymmetric 
state (A1-). 



The relative population of the state Au-|- (upper right panel in figure [TT]) reveals a small 
oscillation in addition to the adiabatic evolution which can be explained by a linear 
stability analysis based on the Bogoliubov-de-Gennes equations (see section [3^ . 

If the time-evolution is initiated with the almost antisymmetric allochtonous state 
(A1-, figure [9]) the wavefunction starts to oscillate between the wells and hence does 
not exhibit a stationary decay behavior so that it is termed dynamically unstable. Due 
to the asymmetry of the system there is no complete population transfer between the 
wells. Because of the infiuence of interaction the Josephson oscillations for small times 
are not sine-shaped but can be described by a Jacobi elliptic function (see e.g. [40] and 
references therein). For longer times the effective nonlinearity decreases due to decay 
so that the oscillation becomes more sine-shaped and its period slowly gets closer to 
the value 27rh/Afig=Q ^ 7.53 for the Josephson oscillations of the linear system where 
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Figure 10. The same as figure IH however, for the allochtonous self-trapping state 
(AIT). 




Figure 11. The same as figure [SJ however, for the autochtonous almost symmetric 
state (Au+). 



A/ig=o ~ 0.834 is the difference between the chemical potentials of the ground and first 
excited state for (7 = (cf. figure [3]). In order to estimate the oscillation period for 
short times we consider the interval between the first two population maxima in the 
left well at ti ~ 4.05 and ^2 ~ 10.84. In an ad hoc approach we take the interaction 
into account by considering the difference between the chemical potentials of the almost 
symmetric, respectively antisymmetric quasistationary eigenstates Au-|- and Al- in the 
middle of the interval [^1,^2] instead of AyU,g=o- In the middle of the interval [^1,^2] we 
still have about 80% of the initial population left so that the effective nonlinearity can 
be estimated as ^2.4 and the respective difference between the chemical potentials 
is approximately given by A/ig=2.4 ~ 0.95 (cf. figure E]). The resulting estimate for the 
period 27r^/A/ig=2.4 ~ 6.6 roughly agrees with the value ^2 — ~ 6.79 observed in the 
dynamical calculation. 

The time-evolution of the allochtonous self-trapping state (AIT, figure [TOj) follows 
the stationary decay behavior until the bifurcation point is reached. After that it tunnels 
completely from the right to the left well and starts oscillating between total population 
of the left well and an intermediate population of the two wells. This asymmetry 
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indicates that a part of the total population of the double well does not take part 
in the oscillation but is trapped in the AuT state. 

In conclusion, we observe that there are dynamically stable states whose dynamics 
is well described by stationary states in an adiabatic approximation. In particular, 
the self-trapping effect is not immediately destroyed by decay but survives until the 
bifurcation point is reached. Especially in the autochtonous self-trapping state (AuT) 
the trapping can be preserved a long time since this state decays very slowly. 



3.3. Finite basis approximation 



Instead of solving the NLSE as a differential equation, one can, in a Galerkin-type 
ansatz, expand the wavefunction ip{x,t) = exp(— i(/i — iT/2)t/h)u{x,t) as 

riB 

u{x,t) = ^Cj{t)uj{x) (28) 
i=i 

using the first tt-b eigenf unctions {uj} and respective eigenvalues {fij — irj/2} of the 
corresponding linear {g = 0) system with the Hamiltonian 



using Siegert boundary conditions. The NLSE can now be written as 

ihdtu = (^Hq — (yU — ir/2) j u + g\u\'^u . 

Since Hq is not hermitian the states {uj} are not orthogonal. Instead they form a finite 
basis set 



(29) 
(30) 



v*Ax) ui{x)dx = 6ji 



(31) 



together with the eigenvectors {vj} of the adjoint Hamiltonian satisfying 

Hlvj = {fi,-iT,/2yv,, (32) 

the so-called the left eigenvectors of Hq (see e.g. [46]). We compute the left and right 
eigenvectors using exterior complex scaling in order to make the states {uj} and {vj} 
square integrable in < x < oo (see Appendix A). 

To calculate the stationary states we set dtUj = 0, 1 < j < tt-b, and consider the 
projections of equation fl5U]) on the left eigenvectors: 

2 rtB 

" , 1 < J < riB (33) 



POO 


riB 


/ dxvUx) 


y^^CjUijx) 


Jo 


i=l 



rtB 

^ClUl{x) 
1=1 

with E = fi — iT /2 and Ej = fij — iTj/2 . Together with the normalization condition 



dx|M(x)| 

I 

and a condition 

arg(M(a)) = 



(34) 
(35) 
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which fixes the phase of the wavefunction we obtain a system of nonhnear equations for 
the coefficients {c^}, the chemical potential /i and the decay rate F. This system fl33|) - 
(135|) is solved with a Newton algorithm. 

Alternatively, equation (133|) can be rewritten as 

"B 



Ci {Ej -E)+^ wi\ c* Cfc Q = 



with 



11 



iki=l 



g I dxv*{x)u*{x)uk{x)ue{x) = wl\ 
'o 



(36) 



(37) 



where the integrations can be carried out once in the beginning. This has some 
advantages for small ub- For the most simple case of only two states this yields the 
nonlinear, nonhermitian 2x2 eigenvalue equation 









with 







Cl 


' + '^wll 


C2 


2 


1 12 * 


f^l2 




C2 


' + 2w\l 


Cl 


2 


1 22 * 
+ W^llCiC2 




= < 


Cl 


' + '^wll 


C2 


2 


1 21 * 
+ W11C2C1 


f^22 


22 

= W22 


C2 


' + Ml 


Cl 


2 


1 21 * 
+ W^22ClC2 



(3^ 



(39) 
(40) 
(41) 
(42) 



In figure[T2]we compare the bifurcation diagram calculated in section [3TT] (solid black 
lines) with the Galerkin approach with tt-b = 2, (dashed red lines) and riB = 30 modes 
(dashed dotted blue lines). Both real (left panel) and imaginary (right panel) parts 




g 




Figure 12. Bifurcation diagram for A(, = 10, \a = 20, b = 1, a = 2. Solid black: 
analytical result (cf. figure [21), dashed red: Galerkin approximation with ne = 2 
modes, dashed dotted blue: ub = 30 modes. 



of the eigenvalues are reasonably well reproduced by such a two mode approximation. 
Naturally the agreement is best for the two lowest states. The Galerkin approximation 
with TT-B = 30 modes almost coincides with the results from section 13.11 in the displayed 
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parameter range < < 3 except for the decay coefficient of the allochtonous self- 
trapping state for high values of g. 




g g 

Figure 13. Bifurcation diagram for = 10, = 20, 6 = 1, a = 2 and projection 
onto tire right eigenvectors {uj}. Solid black: analytical result (cf. figure [2]), dashed 
red: Galerkin approximation with tib = 2 modes, red dots: tt-b = 2 modes and Siegert 
formula, dashed dotted blue: tt-b = 10 modes, blue dots: = 10 modes and Siegert 
formula 

In order to achieve such a good agreement, it is essential to use the biorthogonal 
basis. If only the right eigenvectors are used, i.e. if we project equation ( l30l) on the 
{uj} instead of the {wj}, the results are worse. Figure [T3] reveals that for ne = 2 modes 
the real parts of the eigenvalues are still in good agreement with the analytical results 
whereas the imaginary parts of the eigenvalues clearly differ from the previous results 
since the growth of the decay rates with increasing nonlinearity g (cf. the discussion in 
section [3?T1) is not reproduced. The results for u-q = 10 demonstrate that the agreement 
improves if higher modes are taken into account. Thus the decay rates are much more 
sensitive to the excitation of higher modes than the real parts of the eigenvalues. The 
agreement can be improved if we calculate the decay rates with the Siegert formula ( l26l) 
using the wavefunctions and real parts of the respective eigenvalues calculated with the 
ne = 2 (red dots) and the ub = 10 (blue dots) Galerkin approximation. 

These results indicate that the biorthogonal basis {vj}, {uj} is better suited to 
describe the system than {uj} alone. Naturally the differences between different choices 
of basis sets disappear in the limit u-q oo. 

3-4- Linear stability analysis 

In this subsection we analyze the stability of the adiabatic time evolution of the 
quasistationary states by solving the Bogoliubov-de-Gennes (BdG) equations and 
compare its predictions with the results of the dynamical calculations from section 13.21 
The BdG equations are obtained by linearizing the GPE in the vicinity of a background 
solution ipQ, i.e. by inserting ip = ipo + Sip into the GPE and expanding the resulting 
equation up to ffist order in 6ip. Since we are considering quasistationary rather than 
stationary background solutions we follow the prescription of Castin and Dum [47] 
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for time-dependent background solutions. They argue that only the components Sifj^ 
perpendicular to the time-dependent background solution are relevant. Assuming a 
time-dependence 5il){t) = Sip exp[—{fi — iV /2)t/h — iujt] the BdG equations read in our 



case 



fVjJ 




-9Q*i'f 

-ir/2, where Q = 1 



E 



-H*^P - gQ^li'ol'Q* + E* 




(43) 



with E = /i — ir/2, where Q = 1 — \ipo) {^o\ is the projector orthogonal to the background 
state ipo and the Gross-Pit aevskii Hamiltonian Hqp = Hq + g\'ipo\'^- The population of 
eigenmodes with positive imaginary part of the eigenvalue grows exponentially in time, 
thus leading to instability. The eigenvalue equation (H3ll is solved by expansion in the 
eigenbasis of Hq as described in section 13.31 Figure [TH shows the Bogoliubov excitation 
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Figure 14. Bogoliubov excitation spectra for the same parameters as in figure [4] 
computed with 71b = 4 basis functions. Upper left: autochtonous self-trapping state 
(AuT), lower left: allochtonous self-trapping state (AIT), upper right: allochtonous 
almost antisymmetric state (A1-), lower right: allochtonous almost symmetric state 
(Au+) 



spectra for the same states as in figure H] computed with = 4 basis functions. Note 
that in the figure only seven instead of eight eigenvalues are observed since the eigenvalue 
zero is doubly degenerate. The results do not change significantly if more basis states 
are taken into account. The excitation spectrum of the state Al- (upper right panel) 
has an excitation with zero real part and a large positive imaginary part which leads to 
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instability as observed in the dynamical calculations (cf. figure [H])- The autochtonous 
states (upper left and lower right panel) only have excitations with negative imaginary 
part so that the adiabatic stability seen in the dynamical calculations (figures IHl and [TTi) 
is confirmed. The first (nonzero) excitation of the Au+ state (lower right panel) has the 
eigenvalue uj ~ ±2.1 — i0.0075. The small imaginary part indicates that the excitation 
is only weakly damped. The period 27r/|Re(u;)| 3 of the excitation agrees with the 
period of the small oscillation in the background state's population imbalance (upper 
right panel of figure [TTj) for short times. In the limit g ^ Q the nonzero eigenvalues of 
the BdG equations (143|) are given by the differences of the eigenenergies of Hq so that 
for longer times the oscillation period slowly approaches the value 2'irh/Afig=o ^ 7.53 
for the Josephson oscillations of the linear system. 

The situation is more involved for the excitation spectrum of the allochtonous 
self-trapping state (AIT) (lower left panel). In the dynamical calculation it seems to 
evolve adiabatically until the bifurcation point is reached (figure [10]) . Yet, its excitation 
spectrum shows eigenvalues with positive imaginary parts. These imaginary parts 
Im(u;) ^ 0.0317 are, however, quite small and the characteristic time scale for the 
growth of the respective excitations can be roughly estimated as l/lm.{uj) ^ 32 which 
is on the same order of magnitude as the time At ^ 30 that it takes for the dynamical 
evolution to reach the bifurcation point (cf. figure [10]). Consequently, the evolution 
of the eigenstate is approximately adiabatic up to the bifurcation point if the initial 
population of the destabilizing excitations is small. 

4. Conclusion 

In this paper we analyzed the resonance states and the decay dynamics of the nonlinear 
Schrodinger equation for a double delta-shell potential. 

By means of an approximation the resonance wavefunctions and eigenvalues were 
calculated analytically. In analogy to the the respective closed system the real parts 
of the eigenvalues in dependence on the nonlinearity g show a saddle node bifurcation 
as symmetry breaking solutions emerge for a critical value of g. The imaginary parts 
(decay rates) undergo a similar bifurcation scenario. The approximate analytical results 
are in good agreement with numerically exact calculations based on complex absorbing 
potentials. 

Comparison with a finite basis approximation demonstrates that both left and right 
eigenvectors of the linear [g = 0) system must be taken into account to obtain correct 
values for the decay coefficients with a small number of (e.g. two) basis functions. 

A time-propagation of the eigenstates reveals that the states with linear counterpart 
(autochtonous states AuT and Au+) evolve according to an adiabatic approximation 
based on the quasisitationary resonance states for different values of the effective 
interaction g \ip{x)\'^dx. The time-evolution of the allochtonous self-trapping state 
(AIT) can also be described by the adiabatic approximation until the bifurcation point is 
reached and the wavefunction starts to oscillate between the wells. For the allochtonous 
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almost antisymmetric state (A1-) the wavefunction starts to oscillate immediately so 
that this state never evolves adiabatically. These results indicate that the self-trapping 
is not immediately destroyed by decay but is preserved on a time-scale determined by the 
decay coefficients of the respective self-trapping states. Furthermore it was shown that 
these adiabatic stability properties can also be deduced from a linear stability analysis 
based on the Bogoliubov-de-Gennes equations. 
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Appendix A. Resonance solutions of the linear Schrodinger equation with 
an open double well potential 

The Galerkin-type approach in section 13.31 requires the computation of the resonance 
eigenfunctions u{x) and corresponding eigenvalues /i — ir/2 of the Hamiltonian Hq given 
in equation ( 129|) with the potential V{x) given in equation ( fT4|) which are obtained by 
solving the stationary Schrodinger equation 

-^dl + V{x) ) u{x) = (/i - \T/2)u{x) (A.l) 

with Siegert boundary conditions. The ansatz 

sin(A;x) < x < 6 

u{x) = \ {sin{kb) / sin{kb + d)) s\n.{kx + {}) b<x<a (A. 2) 
{sm{kb) I sin{kb + sm{ka + •&) e^''^^"") x > a 

with k = ^2m(/i — ir/2)/h already makes the wave function continuous sd. x = a,b and 
satisfies the Siegert boundary condition \imx-,oou'{x) = iku{x) as well as the boundary 
condition u{0) = 0. The matching conditions for the derivatives at x = a,b read 

k cos{ka + ^) = {ik - 2Aa) sin(A;a + ^) (A.3) 

and 

k cos( kb) = k ^-^^[^^^ cos( kb + ^) - 2Xb sm( kb) , (A.4) 
sm(A:o + v) 

respectively. The complex quantities k and 'd are obtained by solving these equations 
numerically. The wave function u{x) diverges for x — > oo since Im(/c) < 0. Therefore we 
use exterior complex scaling (see e.g. [1,48]) to make the wave function square integrable. 
The X coordinate is rotated by an angle 6c from the point where the potential V{x) 
becomes zero. In our case this reads 

, . ~ ■ A.5 

" [x — a) exp^ib'cj x> a 



Resonance solutions of the NLSE in an open double-well potential 



20 



In the scaled region the Schrodinger equation becomes exp(2i6'c)M"(x) + /c^M(a;) = 0. The 
equations (IA.3P and flA.4p remain unaltered. For a sufficiently large rotation angle 9^ 
the wavefunction u{x) becomes square integrable in < x < cx). 

Because of Hq = complex conjugation of equation (!32|) shows that the 
corresponding left eigenfunctions v{x) are given by v{x) = {u{x))*. We normalize the 
eigenstates such that dx v* {x)u{x) = 1. 
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